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SUMMARY 


In this contract (NAS8-25101) , a systematic theoretical investi- 
gation of the dynamical behavior of the solar active region has been 
performed. As a result of these studies, we have concluded that the most 
appropriate physical mechanism in helping to understand the disturbed 
solar atmosphere is the propagation of shock waves in a model solar 
atmosphere. During the course of this contract, we have examined two 
important cases: 

(i) The Downward Propagation and Response of the Chromosphere 
In this study, we have examined the responses of the solar 
chromosphere to an infalling material stream resulting from the 
"disparition brusque" of a prominence. We found that the solar 
chromosphere is heated by the shock resulting from the infalling 
material stream and radiation is enhanced. The enhanced radiation 
terminates the shock around the height of the temperature minimum 
in the Harvard-Smithsonian Reference Atmosphere model. This 
radiation enhancement is identified as Optical (H a ) flares. The 
detail of this study was submitted to the National Aeronautics and 
Space Administration, Marshall Space Flight Center, as an Interim 
Report for this contract, dated March 1972 (UARI Research Report No. 
114). A part of this study is also published in Solar Physics , 

Vol. 30, Page 111-120, 1973. 
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(ii) The Upward Propagation of Solar Disturbance and Its Responses 
After completion of the study of the downward propagation of 
a shock through the chromosphere (from ~ 2000 Km to the sun s surface) , 
we felt it logical to examine the responses of the solar atmosphere 
due to an outward propagation shock. Therefore, in this final research 
report, we shall report the results of this study, since the other 
results have been documented already. In this study, we have 
employed the Lax-Wendroff method to solve the set of non-linear 
partial differential equations, because the method of characteristics 
used to analyze the downward propagating shock became invalid due 
to non-homogeneity in the model of the solar atmosphere. It was 
found that this theoretical model can be used to explain the solar 
phenomena of surge and spray. A criterion to discriminate the 

surge and spray was found. The detailed information concerning 

/ 

’ the density, velocity, and temperature distribution with respect to 
the height and time is presented. The complete computer program 
is also included in this report. 

Finally, we would like to summarize the publications and research 
reports resulting from this contract as follows: 
i. Refereed Publications 

(1) "A Kinematic Model of a Solar Flare," Solar Physics , 

Vol. 30, 111-120, 1973. 

(2) "Non-Equilibrium Ionized Blast Wave," J. of Physical 
Soc. Japan , Vol. 36, No. 1, 1974. 

iv 
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(3) 


"Kinetic Description of Solar Wind Interaction with 
’Small 1 Celestial Objects," Rarefied Gas Dynamics 
(ed. K. Karamcheti) , Academic Press, New York, 1974. 

(4) "A Kinematic Model of Surge and Spray," to appear in 
Solar Physics, 1974. 
ii. Research Report 

"Propagation of Downward Shock Waves Generated by Infalling 
Dense Prominence Materials in a Realistic Solar Atmosphere," 
UARI Research Report No. 114, March 1972. 


v 



TABLE OF CONTENTS 


CHAPTER 

ACKNOWLEDGEMENT 

SUMMARY 

I INTRODUCTION 

II FORMULATION OF THE PROBLEM 
II-l Hydrodynamics Model 

II- 2 Initial and Boundary Conditions 

II-2-1 Initial Conditions 

II- 2-2 Boundary Conditions 

III NUMERICAL TECHNIQUE 

III- l Conservations! Form of Equations 
III-2 Numerical Stability 

III- 2-1 Time Interval 
III-2-2 Artificial Viscosity 
III-2-3 Shock Dissipation 

III-2-4 Minor Modifications on Difference 
Equations 

III-3 Computation Procedure 
IV RESULTS OF COMPUTATION 
V CONCLUSIONS AND RECOMMENDATIONS 
REFERENCES 
LIST OF FIGURES 
APPENDICES 


PAGE 

ii 

iii 

1 


4 

4 

7 

7 

9 


11 

12 

16 

16 

17 

19 


21 

24 


26 


30 



CHAPTER I 
INTRODUCTION 


During previous studies [I, 2] » we found that the optical (H^) 
flare can be identified with the response of the solar chromosphere to 
a shock wave propagating downward through the chromosphere. The shock 
wave is related to an infalling material stream resulting from the 
"disparition brusque” of a prominepce as suggested by Hyder [3, 4, 5], 
and Nakagawa and Hyder [6]. In the general impact theory, there are 
always two waves generated propagating in opposite directions right at 
the moment of the impact. According to the coordinate system which 
we adopted here, one of these two waves is propagating downward through 
the chromosphere to the photosphere, and another one is propagating upward 
through the transition region to the corona and beyond. The study of a 
downward propagating shock through the chromosphere has been, completed and 
reported [1, 2]. Therefore, we shall present the results of upward 
propagating disturbances (either shocks or subsonic disturbances) 
through a model solar atmosphere in this report. 

To calculate the downward propagating shock through the chromo- 
sphere, we have used the CCW (Chisnell, Chester, Whitham) [7] approxi- 
mation which is based on the theory of characteristics. However, it 
is noted by Bird [8] that this method of approximation has only provided 
satisfactory accuracy for a shock propagating into a denser medium. 
Therefore, it seems to us that it is not quite proper to apply the CCW 
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approximation to calculate the outward propagating disturbances, because 
the model so-lar atmosphere is attenuated outward. Thus, we have chosen 
the Lax-Wendroff method [9, 10] for the present study. This is a 
numerical method, which has the advantage of taking care of both the 
sub- and supersonic disturbances. The numerical accuracy can be con- 
trolled by using proper numerical techniques, such as by specifying the 
proper time increament and grid size in the computation processes based 
on the physical model of the solar atmosphere. 

In this study, the evolution of the disturbances, originating at 

5 

~ 30,000 Km (*“ 0.0428 R s R s being the solar radius ~ 7 x 10 Km), 
is examined in detail using the method we mentioned in the previous 
paragraph. These disturbances are identified as pressure pulses with 
different strengths and durations. From the present results, we have shown 
that the short duration (few minutes) and moderate strength pressure 
pulse (A p ~ 2) will result in the phenomenon of "surge," because it 
shows that a stream velocity of -100 - 200 km/sec can be achieved in this 
case”! It also demonstrated LliuL tliuie lu an euoeutd*aX- 4 Lax£_rif . the material 
falling back to the surfs surface which agrees with the observation. 

The longer duration (i.e., 20-40 minutes) and stronger strength pressure 
pulse (Ap ~ 10) shows that a stream velocity of the order of 1000 km/sec 
is achieved and no falling materials can be seen. Thus, we have identified 
this case as the "spray". 

Finally, we should point out, that there is little difference between 
the adiabatic calculation and the calculation with a Cox-Tucker type 
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radiation loss as has been shown in the calculation of the downward 

propagating case [1], This is because the Cox-Tucker [11] radiative 

loss is mainly based on the hydrogen-equilibrium radiative-equilibrium 

estimation that only covers radiation such as the Balmer 

series. It is worthwhile to examine the radiative effects in further 

detail. 
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CHAPTER II 


FORMULATION OF THE PROBLEM 
II-l Hydrodynamics Model 

In the solar atmosphere of interest in the present study, the 
gyro-radius is ~ 10 km, and while the scale height is of the order of 
a thousand kilometers, we can consider that the medium is filled with 
collision-dominated plasma. Thus, the physical behavior of this plasma 
can be considered as a continuum fluid. Consequently, the hydrodynamic 
model was chosen for the present problem. 

In dealing with radiative cooling effects in this problem, we 
have chosen the Cox-Tucker model [11], because the dominating radiation 
in this part of the solar atmosphere results from brerasstreahlung , 
recombination radiation, and collision- induced line emission. A summary 
of this radiative loss is shown in Figure 1. For the convenience of the 
numerical calculation, a simple, analytical closed-form expression is 
adopted, 

Qr = X P 2 T a • (2-1) 

The symbol is the radiative cooling rate (ergs /cm 3 sec), p is the gas' 
mass density (gm/cm 3 ) and T is the temperature CK). Finally x and a, 
constants determined by the results given by Cox and Tucker, as shown 
in Figure 1. The numerical values for these two constants at various 
temperatures are presented in Table 1. 
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Table I. Temperature of Radiative Cooling Rate 


Range of Temperature °K 


T < 5.0 x 10 4 


5,0 x 10 4 < T <. 2 . 5 x 10 s 


2 , 5 x 10 s < T < 7. Ox 10 6 


7 x 10 6 < T 




a 

a - 

3.55 

a = 

0 

a = 

-1.172 

a. = 

0.288 

T a (ergs/ 

cm 3 sec) 


X=* 3.0 x 10‘ 


x = 1.0 e 3Z 


X = 1.0 X 10 : 


As we have discussed previously, the hydrodynamic model can be 
used for the present study. In order to avoid unnecessary complexity 
while retaining the basic physical process of the problem, we consider 
the plasma flow guided upward along a vertical magnetic flux tube. The 
possibility of such a confinement is discussed by Nakagawa and Hyder [12], 
and it was shown that confinement is possible when the gas pressure 

within the plasma flow is smaller than the local magnetic pressure, i.e.. 



— > JX. m 2 p , 

8tt y+1 


( 2 - 2 ) 


where B is the strength of magnetic field induction, y the ratio of 
specific heats, M the shock Mach number and p the gas pressure. If we 
consider that the maximum gas pressure in the model of the solar atmo- 
sphere is - 1. 5 x 10“ 1 cgs , which corresponds to the gas flow velocity 

of the order of 4 x 10 2 Km/sec , we find that Eq. (2-2) is satisfied for 
B 60G, which is a reasonable value of B in an active region. Therefore, 
the magnetic force effect can be ignored in this calculation. The 
governing equations for the present problem can be written as 
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Continuity : 


t 


7“ + V • (p V) = 0 , (2-3) 

0 t 

Momentum: 

p fj + V • V\? = - Vp + p| , (2-4) 

Energy: 

^ <\ c -*■ -> 

p + (v V) e = -pV*V + V(KVT) - Q r (2-5) 

3t 

where V is the flow velocity and g is the gravitational acceleration along 
the normal axis from the sun's surface, toward the sun. 0 is the internal 
energy of the gas per unit volume, K is the thermal conductivity and Q R 
is the radiative cooling loss rate given by Eq. (2-1). Finally the 
equation of state is 



( 2 - 6 ) 

(2-7) 


with R, k, m, and C v being the gas constant, Boltzmann constant, average 
molecular weight and specific heat at constant volume, respectively. 

Let us adopt the spherical coordinates for the present study, and 
further assume that the case of spherical symmetry, the thermal con- 
duction is negligible compared with radiation. The Eqs. (2-3) through 
(2-5) become 

!£ + ■£- <Pu) + — = 0 (2-8) 

dt dr r 
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3u 3u 
Tt + U 3r 


IiE - 


(2-9) 


3e 3 e 

37 +u S7 


£ te + uil 

p \3t 3r. 


where u represents the radial velocity and depends on radius r 


and time t - 


The gravitational acceleration g is given by 


Ss Rs 


with g being the gravitational acceleration at the surface of the sun, 

D 

and R„ being the solar radius. 


Eq. (2-8) through (2-10) 


are set of non-linear time dependent 


partial differential equations without dispersion coefficients (i.e., 
viscosity, diffusivity, etc.). To find an analytical solution for this 
set of equations with various boundary conditions is impossible. How- 
ever, it is possible to obtain a numerical solution. Some discussion 
on the existing numerical method will be given in the next chapter. 

II- 2 Initial and Boundary Conditions 
II-2-1 Initial Conditions 


Initially, we have assumed the solar atmosphere at the photosphere 
and chromosphere to be the Harvard -Smithsonian Reference Atmosphere (HSRA) 
[13], Beyond its range, the solar atmosphere is assumed to be in a state 
of hydrostatic equilibrium, which can be calculated from Eq. (2-4), thus 


« - p° g = _ p° g g Rg/r 2 , 


(2-12) 
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where superscript o denotes the quantities at steady state. 


Substituting p° = p°RT° into (2-12), we obtained 


dp°(r) 

dr 


- P°(r) 
o 


g S R S 

RT°(r) 



1 

T° (r) 


dT°(r) 

dr 


(2-13) 


Eq. (2-13) should be numerically integrated for a given temperature 


profile of the atmosphere. However for an isothermal atmosphere Eq. (2-13) 
reduces to a simpler form 


1 dp° __ _ &s R s 

p 0 dr R T ° r 2 


(2-14) 


Integration of this equation yields the solution for the density profile of 
the hydrostatic atmosphere, 


p° = p° Exp 
o 


R s (l 1 V 

L RT ° \ r ' r i/. 


(2-15) 


where p is the reference density at r 

o 

distribution is then 


r l* 


Steady state pressure 


o o 
P = P Exp 
o 


r§ r 2 

B s s 

/i - 

JL\" 

RT° 


r J. 


(2-16) 


where p° is the pressure at r = r^. 
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II-2-2 Boundary Conditions 

There are two boundary conditions, one at r - 0 and one at r 
to be given, The lower boundary condition will be characterized by 
disturbances such as density pulse, temperature pulse, velocity pulse, 
and pressure pulse. To introduce these pulses, we can specify them by 
prescribing the amplitude and duration of the pulse which depends on the 
characteristics of the disturbances. 

For example, the disturbance is introduced at the lower boundary 
as a pressure pulse and the velocity is zero at t = 0. But, the velocity 
on the boundary for t > 0 will be determined from the continuity equation 
thus 

u n » p n u n /p n (2-17) 

1 2 2 1 

where subscript 1 denotes the boundary and 2 denotes the mesh point 
next to the boundary, superscript n denotes the time increment, such 
that t “ nAt < T, T being the duration of the disturbance. When t > T, 
disturbance is gone, and the lower boundary returns to its unperturbed con- 
dition, i.e., a hydrostatic equilibrium state. 

On the upper boundary we have used the completely absorbed con- 
dition, i.e., all the effects due to reflected waves are ignored. The 
reason for keeping this assumption is because the time for the reflected 
wave to reach the upward-propagating wave packet is much longer than 
the time for the upward-propagating wave to reach the upper boundary. Thus, 
the non-ref lective boundary condition for velocity can be expressed as 
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where $ represents a physical quantity, the subscript j denotes the raefih 
points at upper boundary and superscript n denotes the time step. 



CHAPTER III 
NUMERICAL TECHNIQUE 

Numerical computation of time dependent inviscid, compressible 
flow is a formidable task because of the appearance of discontinuities 
in the flow field. The conventional way to solve hyperbolic types of 
differential equations is the method of characteristics. In this method 
partial differential equations are usually written in characteristic form. 

Due to the presence of discontinuities in the fluid, these characteristic 
equations can not be integrated over the entire region of space. Instead, 
the integral form of the differential equations is used for the dis- 
continuities while differential equations are applied to the remaining , 
region. Although this method is very simple in principle, its application 
to practical problems is very lengthy and cumbersome. Furthermore, the 
fact that one cannot know the time and location of such discontinuities 
in flow field prior to the computation makes application to actual problems 
almost impractical. Watts and Rosenverg, et.al. [14] solved the transient 
adiabatic compressible fluid flow in a duct by using the characteristic method 
in an elegant manner. Their method can be, in principle, extended to 
solve the present problem. However, the resulting computational procedures 
will be too complicated to be practical. It may be worthwhile to examine 
this method in some detail at a later date. Gentry, et al., [15] 

used the FLIC method, known as Fluid in Cell, to describe the time dependent 
equations of motion for the compressible flow of a fluid. This method has 
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been used to solve a wide variety of problems in compressible fluid flow. 
Hundhausen, et al, [16] used this method to simulate the flare generated 
disturbances in the solar wind. 

Currently, the most commonly adopted method for solving the 
compressible fluid flow problem is probably the Lax-Wendroff difference 
method. Lax-and Wendroff [9, 10] suggested that partial differential 
equations are first written in divergence free form, and then the 
difference equations in the divergence free form can be generated from 
these equations. The mathematical proofs are beyond the present scope 
of these studies and will not be presented in this report. However, the 
basic idea of this method is that errors caused by the discretization 
process tend to smooth the solution. This allows the representation 
of shocks by smearing discontinuities over several mesh points. 

There are many versions of the difference scheme for the conser- 
vational form of equations. Some detailed comparisons among these versions 
are made by Ehmery [17] and Burstein [18, 19]. One version due to Burstein 
[20], is used for the present study. 

III-l Conservational Form ■ of Equations 

Divergence-free form of the governing equations will not create or 
eliminate any flow variables. After some algebraic manipulation (see 
Appendix A) the governing equations [2-8], [2-9], and [2-10] can be written 
in Eulerian pseudo-conservational form; 
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(3-1) 



(3-2) 


(3-3) 


where E is the total energy per unit volume, given by 


pu‘ 

E - + T 


(3-4) 


Using vector notation, Eqs. (3-1), (3-2) and (3-3) can be put in the 
form 


ML - -M +K 

3t 3r 


where U-, F and K are three components vectors; 



(3-5) 


(3-6) 



(3-7) 
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(3-8) 


2pu 

r 


K = 


-Pg 


_ 2£U : 


-pug - Q B - 


j^uCyE - (Y-l) ^5“)] 


Using the difference operator. Equation (3-5) is approximated by 


6 t U = -6 F + K, where (3-9) 

6 t and <5 r is yet to be discussed. There are several versions of Lax-Wendrof f 
difference schemes which have been extensively used for a wide range of 
fluid flow problems. 

Lax-Wendrof f scheme is based on the Taylor series expansion of the 
vector function U(r, t + At) so as to include the second order term 
3 2 U/9t 2 . A two-step method obtained by Richtmyer [21] is used here. 

The values at the intermediate points are computed at a time t + At/2 
using a first-order correct scheme, and then a second-order correct scheme 
(leap frog) is used to compute the value at time t + At. The overall 
scheme has, then a second-order correct differencing scheme. Burstein, e t « j a l. 
[20] suggested that instead of computing the intermediate value at time 
t + At/2, they compute them at t + At and then average the F difference 
at t and t + At so that both the U and IF values are centered at point 
(i, t + At/2). The above differencing method applied to Equation (3-5) 
yields the following difference approximation: 


14 



Intermediate Values 


U n+1 = i 
~i+l/2 2 


U n + U n 1 - 
'i+l ~i 





— n+1 
~i-l/2 


= ± / un 


uy + uy 

~i-l 


At 


_n 
F . - 


? n 

' 1-1 


+ 4Mk" + k” . 

2 V^i ~i-l 


(3-10) 


U, 


’n+1 


T n 


T n 


Vi + Vi 


At 


2Ax 


.n 


,n 


F” - F 
“"I+l ~i-l 


+ At K' 


n 


where the bar signifies the intermediate flow variables at (i-l/2)Ax, iAx, 

(i+l/2)Ax and at (n+l)At. Using these intermediate values of the variables 

we calculate the final values. 

Final Values 


n+i 

~i 


U?" 

~i 

At 

1 

- 

l~i+l 

2 Ax 



* ( 

n 

K 

~i 

-n+1 


7 n 
-i-1 


— n+1 

where the intermediate values of £,±+ 1/2 


/—n+1 
+ (~i+l/2 


- Cm) ] 


(3-11) 


— n+1 — n+1 

ii- 1/2 and are 


calculated by using Eqs. (3-7) and (3-8) with intermediate flow 
variables obtained in the previous step. Combining Eqs. (3-10) and 


(3-11), it is easily seen that this is a second-order correct -differencing 
scheme . 
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III-2 Numerical Stability 


Finite difference equations may exhibit rapidly growing and 
oscillatory solution that cannot resemble the true solutions of partial 
differential equations. In this case, the difference equation is said 
to be computationally unstable. The origin of instabilities varies with a 
particular set of partial differential equations. Many theories and 
criteria have been developed by many investigators, such as Richtmyer 
and Morton [22], Hirt [23], and Van Leer [24]. However, . it is not possible 
to have a general theory for higher-order non-linear equations, such as in 
the present problem. The discussion given here is not rigorous, but 
presents some heuristic techniques which prove very useful in stabilizing 
the computations. 

III-2-1 Time Interval , 

As a first approximation, the time interval for successive interations can 
be found by applying the stability criterion of Cour ant, Friedrichs and Lewy, 


At < 
s 


Ar 

u | + a 


(3-12) 


where a is the local sound speed. When a thermal conduct ion term is included 
in the governing equations, thermal conduction stability criteria should 
be considered. Applying the Fourier method proposed by von Neumann con- 


cerning the heat conduction equation, one finds the time interval [5j At c to be 


At c < 


j P (Ar) 2 
KCy-DT 7 / 2 


(3-13) 
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In the case of heat exchange terms becoming arbitrary, namely Qr , the K3 
vector in Eq. (3-8) of the present problem, can be varied very rapidly. . 
When such a case occurs, numerical oscillation starts and computation 
terminates. In order to smear out such oscillation, the time step is 
chosen according to 

I Q r! At Q < t • < 3 - 14 > 

For a model that does not include the conduction term, 0.25 At is 
sufficiently small to satisfy all stability criteria. 

111-2-2 Artificial Viscosity 

If the strength of the disturbance is large, a sharp jump in the 
flow variable occurs in the flow field. We can find this jump condition 
by applying the Rankine-Hugoniot relations. In numerical simulation, these 
discontinuities easily cause the onset of instability. Thus, the artificial 
viscosity is introduced to help offset the instability. The idea of 
introducing artificial viscosity into shock calculations is due to the 
work of von Neumann and Richtmyer [25]. The basic requirements for a 
purely artificial dissipative term are: 

(i) all flow variables should have smooth transitions across the 
discontinuity; 

(ii) transitions should have correct jump conditions computable 
with Rankine-Hugoniot conditions; 
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(iii) the discontinuity should travel at very nearly the correct 
speed ; and 

(iv) the thickness of transition is independent of shock strength, 

pressure or density of material while the shock is moving [26], 
Introducing a pseudo-viscous pressure term in the compression zone, 
it is shown that all requirements for the artificial viscosity are 

t 

satisfied. The pseudo-pressure is given by [26] 


q « 


2 

(pi 2 ) if W3r < 0 

0 if 3u/3r > 0 


(3-15) 


where A is a constant having the dimensions of length. Then the total 
energy, Eq. (3-4), for the compression region, 3u/3r < G, is modified 
to include the pseudo-pressure, such that 


E = P-+ <L + PHI 

Y - 1 2 


(3-16) 


It is seen clearly that this correction only affects the compression region 
and that the continuity equation is inact by this modification. Letting 
i - $Ar s the appropriate difference approximation for the altered pressure 
is then 




P? 


l 




for the intermediate step, and 


if 


n n . 
u i+i <u ± 


if u 


n 

i+1 



(3-17) 


(p+q) 


n+1 
i+1/ 2 


P 


n+1 + i 
i 2 


P 


n+1 

i 


3 2 


.n+1 

1 + 1/2 


n+1 n+l\ n+1 n+1 

u i+l - u 1+ l/2 / lf 5 i+1 < S i+l/2 

_ n+1 _ rr+1 

lfu l+l — U i+1/ 2 
(3-18) 
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for the final step of computations, where — signifies the intermediate 
values, and 8 has the value of 1 to 5. The effect of q on the overall 
picture is carefully tested by several trial runs. There is no significant 
change in flow variables except at discontinuities, and the transition 
occurred over 3 or 4 mesh points. 

III-2-3 Shock Dissipation 

Mechanical energy carried by the shock wave is dissipated into thermal 
energy of the gas which experiences an irreversible, non-isentropic process 
as the shock front passes through the gas. For a unit mass of gas, the 

thermal energy increases in terms of enthalpy which depends on how the 
post shock gas returns to its pre-shock gas state. Schtzmann [27] suggested 
that the gas expands adiabatically until it comes back to the initial 
pressure and then cools down until it reaches the initial density. Along 
this path, the change of enthalpy is given by [28] 



where subscript o denotes the pre-shocked gas. Then, for a periodic disturbance 

propagation, the total enthalpy change of the gas is 

Q d = - Ah (3-20) 

w 

where u) is the period of the disturbance. Letting 

= £_, a = £- and T = (4) 1/2 , 

Po p o o 

Eq. (3-20) becomes 
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(3-21) 



(o y 




x ) 


a 2 +D 

Y-l 


* 


where \p and 0 are functions of shock Mach number; that is, 


=. 21- M 2 - III 

Y+l s Y+l 

(Y+l) M s 2 
(T-l) M g 2 + 2 


(3-22) 


shock speed relative to pre - shock gas 
s sonic speed in pre -shock gas 


(3-23) 


The detailed structure of the shock front must be known in terms of its 
position and its strength M g in order to find the accurate shock dissipation 
in the gas. However, it is not possible to determine the exact position 

of the shock front, since the transition occurs over several mesh points. 

It is not quite clear how to determine the exact shock strength M g in 
this numerical calculation. The pressure difference between the neighboring 
points does not give the shock strength, because the unperturbed solar 
atmosphere possesses a densiyt gradient. In order to give an 

approximate shock strength at position iAr, the following equation is 
employed: 


M 


s 


i 


u . i — u , 

i-l i 


a , 


(3-24) 


i 

This equation only gives a parameter which is related more to the local 
gas flow than the shock strength. In other words, Eq . (3-24) is 
a sufficient condition for a shock, but not a necessary one. It does. 
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however, provide a mechanism for continuously monitoring the presence of the 
discontinuity, whenever the discontinuity occurs, in a very simple manner. 

In computation, then, from Eq . (3-21) with Eqs. (3-22) and 

(3-24) , the dissipated energy due to the shock is included in the form 

r 1 



Pi a i 
0) y 


0 





T 

y-1 




for £.1 (3-25) 
i 


For a dynamic model, the exact value of w cannot be defined. An estimation 
of a) is made on the ground that the weak shock travels one mesh point with- 
in the time At_, i.e., 
s 


At s 


Ar 

u + a 


Since the Lax-Wendroff method is an explicit difference scheme, disturbances 
travel one mesh point for each full iteration. Thus, 0 ) equals to At- 
approximately. Without the addition of Q 0 in the K 3 term in equation (3-8), the 
temperature of the gas just ahead of the shock. goes negative sometimes,' and the 
calculation is terminated. With this modification, the solution remained 
stable and there is no noticeable difference in shock structure. 

I II- 2-4 Minor Modifications on Difference Equations 

Due to the exponential decrease of density in a quiet solar atmosphere, 
the difference scheme needs two minor modifications. For a hydrostatic 
solar atmosphere, the density variation along the r-direction is 
from Eq. (2-15), 


n ~ n . n , . v 

P i = P 1 SXP (“ h i r i/ Ar ), 


(3-26) 
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where 


n 


RT 


n 


Ar (scale height) 


(3-27) 


Then the nearest- neighbor average value appearing in the first term 

in RHS of Equation (3-10) is larger than the value at the central point, i.e., 


0.5 fp n + p n \ =0.5 p 1 ! (e“ h i + e ) 

V i-1 i+1/ L 1 v . 


■ °- 5 P I [ 2+ ( h i) 2 


nV 


p n u+ 

1 V 2 


> p 


n 


(3-28) 


If this difference is not corrected, the density after the intermediate step 
will become excessively large, and the unchanged vertical pressure cannot 
support this excess material. Consequently, a downward velocity appears over 
the entire field. Instead of a simple average, an expression 


0.5 [p n + p n 

i-1 i+1. 


\ + 

p* 1 - 0.5 

/V 1 , 

+ P 1 ? 

/ 

i 

y i-1 

i+iy 


is used for 


Rearranging this, one gets 


0.5 


+ p n ^ + 0.5 p n 
i+1/ i 


1-0 


. 5 (e' h i + e h i) 


(3-29) 


If the nearest-neighbor points have the same value as the center, then 
Eq. (3-29) reduces to the simple neighbor-points average value. 

In the process of calculation, the 2nd correction term in Eq. (3-29), is 

applied for the intermediate step of each iteration for all conservational 

variables at full mesh points. 

Due to the exponential variation of the hydrostatic equilibrium state. 
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the difference approximation is systematically different from the 
derivatives they approximate. That is why the error grows from Eq. 


(3-26), 




(3-30) 


But the simple centered difference approximation to the first derivative 
with respect to r gives 


(3-31) 

K 

dr 

Unless this discrepancy is corrected, an incorrectly calculated pressure term 
in Eq. (3-2) will set the flow field in an upward motion. In order 

to avoid this non-physical situation, a correction term c" is defined such 
that 

( 6 r + C l) P i - ” (3 ~ 32 > 

Combining Equations (3-31) and (3-32), c™ is found to be 

n L 

C i Ar 

This correction term is included for both steps in each iteration. 


h n + 0.5 


e - h ± _ e +h 


*) 


(3-33) 
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III-3 Computation Procedure 

The actual calculation procedure for the present problem is 
illustrated in Figure 2. 

The steady state temperature T of the atmosphere is assumed to be known 
and the steady state density is found by using Eq. (2-15), i.e.. 


n n _ 

p = P Exp 
i 1 


g„ R 


's s l 1 
r i 


RT n 

i 


r l 


(3-34) 


For a hydrostatic equilibrium state the velocity of the field is zero, 


n . 

u = 0 

i 


(3-35) 


The total energy is then, from Eq. (3-16), 


n 


p? + + (p 


y - 1 


2P- 


(3-36) 


n xi 

where p and q are determined by Eqs. (3-17) and 
i i 

(3-18), respectively. 

The disturbance at the lower boundary is introduced in terms of a 

pressure jump. For instance, p*? — 2p**, T n = T n will give p n /p T1 — 2. 
r * 1*1 *1 

This arbitrary disturbance was kept constant for a prescribed time 
interval t. After this period of time, the lower boundary returns to its 
original state. 

The time increment At is found by the CFL (Courant-Fredricks-Levy) 
condition, and is applied to each mesh point in the flow field. 
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Using these initial values of conservational variables, the 


fluxes at time (nAt) are found from Eq. (3-7) and (3—8) 


-n>n _ 

V 


n 


(pu) 


n 


(Y-l) E" 


- ^ »; <«;>■ 


^i 


K 


n 


U n ( y E n 

i \ Y \ 


2P5 Uj 


(I"!) n ,„n 


p” 


K, 


n 


K 


n 

3± 


where 


o" * - 2 P ? (U ? )2 

H i S i 

Y 


,_n n „ , n \ n ±. 
p i u i 8 i ” Vi ” r 


n / n , PS ( u ?> 2 

“i W E r (Y-i) — ; 


+ (Qd) 


g^ = g s R s a /(r^) 2 ; gravitational acceleration 


E n is given by Eq. (3-25) 
i 


<v” - 


- x p" (T")“ 


n 


and (Q d ) is given by Equation (3-25). 

Using these fluxes at nAt, the intermediate conservational flow 
variables at (n-KL)At are found from Eq. (3-10) with Eq. (3-29). Inter- 
mediate fluxes at (n+l)At are found by similar manner as those at nAt 
except intermediate flow variables, U^ 1 should be used in Eq. 

(3-2 9) and (3-33). The new flow variables at (n+l)Atare then calculated 
by using Eq. (3-11) with the aid of Eq. (3-33). 
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CHAPTER IV 


RESULTS OF COMPUTATIONS 

Numerical results are obtained for given various initial boundary 
conditions identified as a pressure pulse with different amplitude- and 
duration resulting from possible solar disturbances due to solar activities. 

All disturbances are placed at the lower boundary which is located at 

-'30,000 Km (1.043 R s ) above the sun 1 s surface, and all the calculations are carried 

out to ~ 3 R s (R s being the solar radius). The results obtained in this 

report are the density, temperature and mass flow velocity as a function 

of height and time for Ap (pressure disturbance) equal to 2 ~ 10 and 

At (duration of the disturbances) equal to 30 sec, 60 sec, 120 sec, 1200 

sec, and in some cases, 2400 sec. A detailed discussion of these results 

will follow. 

Figure 3a, b, and c plotted the disturbed density, temperature and 
velocity due to disturbances of Ap = 2 and At = 120 sec. It shows that 
the disturbance has little effect on temperature and density and its 
influence on mass flow velocity is significant. It appears that this 
disturbance has created a mass stream shooting out from the upper chromosphere 
or/lower corona to upper corona with a flow velocity ~ 50 Km/sec at ~ 1R S 
from the surface (i.e., ~ 2R S from the center of the sun). Similar plots 
with different initial strength and duration of disturbances are 

given in Figures 4a, b, c through 11 a, b, c. These results clearly 
demonstrate that the characteristics of the disturbances are the essential 
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parameters of the effects of the disturbed solar atmosphere. 

Some general features of the disturbed solar atmosphere can be 
observed from these results. Namely, the stronger initial disturbance 
gives a stronger temperature enhancement and its mass flow velocity can 
reach as high as 1000 Km/ sec , and a longer duration of the disturbance 
will sustain the disturbed solar atmosphere, and there will be material 
falling back to the sun's surface. For example, we have plotted the 
velocity versus height for Ap = 6 and At = 30 sec, 120 sec, 1200 sec, 
and 2400 sec at t - 40 min. after explosion. This shows that the negative 
velocities (i.e., downward velocity) appeared near the surface of the 
sun for At. = 30 sec, 120 sec and 1200 sec and for At = 2400 sec, mass 
flow velocity just ejects out from the sun's surface all the way. From 
this evidence, we may suggest that the surge develops due to a short 
duration disturbance, because, observations show the material falling 
back to the sun's surface during a surge. - 

\ 

Now, we shall calculate the total energy of the disturbance initially 
introduced into the solar atmosphere. The total energy can be computed 
from 


E = & 


Y+l 


(E V iPl u?) 


(4-1) 


where on the right hand side, the first term represents the internal energy 
and the second term represents the kinetic energy of the gas in a volume 

(IV.) . The results corresponding to various disturbances are shown in Table II, 

i - ' 


27 



Table II 

Total Energy Per Cross-Section Area for Each Disturbance 

(ergs/Km 2 ) 


Ap 

At = 30 sec 

At - 120 sec 

At - 1,200 sec 

At - 2,400 sec 

10 

7.31 x 10 17 

9.42 x 10 18 

3.80 x lO 20 

6.72 x lO 20 

6 

1.98 x 10 17 

2.63 x 10 ie 

1.38 x lO 20 

2.76 x 10 20 

2 

2.60 x 10 16 

1.48 x 10 17 

1.8 x 10 19 



If we consider the cross-section area of a disturbance which has a radius 
of " 500 Km., it will give a total energy of ~ 5.28 x 10 26 ergs for the 
Ap = 10 and At = 2400 disturbance. This may correspond to the total energy 
of a class of sub-flare. 

From those density profiles, such as Figures 3a, 4a, 5a, 6a, 7a, 

8a, 9a, 10a, and 11a, we can estimate the amount of particles which can be 
ejected into the corona (or solar wind), and the results are given in Table III . 


Table III 

Total Number Particles Per Cross Sectional Area 
for Each Disturbance (///Km 2 ) 


Ap 

At = 30 sec 

At — 120 sec 

At = 1,200 sec 

At = 2,400 sec 

10 

4.77 x 10 26 

5.87 x 10 27 

1,87 x 10 29 

3.44 x 10 29 

6 

2.60 x 10 26 

3.37 x 10 27 

1.42 x 10 29 

2.85 x 10 29 

2 

5.5 x 10 25 

7.30 x 1G 26 

3.43 x 10 2s 



28 





































Again, if we consider the cross-section area of the disturbance being 
~50Q Km in radius, we find that 2.7 x 10 35 particles can be added to the 
corona, which is believed by many to be a reasonable number. 

From this study, we have shown that the surge and spray can 
result from disturbances in the solar atmosphere. After the 
disturbance has been introduced, the corona may settle into a new 
equilibrium state. Evidence for this has been reported in some of 
the observations from the ATM/Skylab experiments. 


29 



CHAPTER V 


CONCLUSION AND RECOMMENDATIONS 

In this investigation, we have examined the upward propagating 
solar disturbances in a model atmosphere. It was found that the 
characteristics of the disturbances have dominant effects on the 
disturbed solar atmosphere. We may conclude from this study that the 
phenomena of surge and spray can be discriminated by the characteristics 
of the initial disturbance, as we discussed in the previous chapter. Also, 
the present model. can be used to examine the observed X-ray data from the 
Skylab mission by relating the X-ray emission to the dynamical responses 
of the solar atmosphere. The initial disturbances introduced in this study 
can be either subsonic or supersonic without limitations. 

The radiation effects on this problem were examined by using the Cox- 
Tucker radiation loss function. We found that there is no noticeable 
difference between the adiabatic calculation and the radiative calculation 
with the Cox- Tucker radiative loss function. This is due to' the fact that the 
radiative loss function given by Cox-Tucker is decreasing as the 
temperature is increasing. Therefore, it is necessary to calculate the 
radiative loss energy from the spectral lines in order to have more 
accurate results. Also, we have ignored the transport effect in the present 
analysis. 

Finally, we shall outline as follows, the steps which should be 
taken to improve the present analysis;. 

(1) Include magnetic field in this model calculation. 


30 



(2) Include thermal conduction effects. 

(3) A detailed radiative hydrodynamic calculation 
procedure needs to be considered. 
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APPENDIX A 


Derivation of Conservational Form of Equations 

The governing equations for the present problem written in 
spherical coordinates are 


Continuity; 


|£ = - 3 (pu) - IE- . 
dt dr r 


(A-l) 


Momentum; 


3u 

3t 


3u _ 3. _3P_ 

3r p 3r ^ 


(A- 2) 


Energy; 


p — = - pu — + ^ + u - Qp, but by virtue 

3t 3r p \3t 9r/ ^ 


of Equation (A-l), this can be written 


3e 3e 3 u 2u 

p 3t = _pu i7‘ p 97' p T ' 


(A-3) 


Continuity equation (A-l) is already in conservational form. 
Multiplying p to equation (A-2) and rearranging with the aid of equation 
(A-l), we have momentum conservational form of equations, 


iM>= _ JL (P u 2 +p) - Pg - , 

dt dr r 


(A- 4) 


noting however p = (y-l)E - pu 2 by definition, equation (A-4) becomes 


3(pu) _ _ ± | (y _ 1)E _ Qt3) pu2 l _ pg . 2Pi 

J r 


3fc 3r 1 " 2 

Energy conservational form is found by adding four equations, i.e., 

, x eqn. (A-2), — x eqn. (A-4) and eqn. (A-3). 

2 2 


(A- 5) 


£ x eqn, (A-l) 



Collecting 3/3t terms, we have after some simplification, 

JL (pe + fid) - |£ (A- 6 ) 

3 1 2 3t 

Collecting 3/3r terms and other remaining terms, we get, after 
some simplifications, 


■ 37 [ u(p + ° E + T -5 ] - “Pg - T [p £ + ^ + p] - \ 

(A-7) 

“ - 3 ^ [u ( TE - - Pu z )J - pgu - Q r - - Ju(rE - pu 2 )J 


Thus, energy conservational form of equation is 
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98* C GS«GRaVJTaTIONaL acc»at The surface OF the sun 

6 3* C DX=SPaTIAL INCREAMENT 

9.6* .. C_ RsG A S CONSTANT ' . 

97* C 0 T # U M E INCREAMENT 

98* C D T P » T I M E INCREAMENT TO RECORD Da^aS ON THE MAGNETIC TAPE 

99* ... .C UlCaUwpFRTURBF.O. CORONA OENSlTr AT THE SuR F Ac £ O.f THE SUN 

S'.!* C V A R I ABE NAMES 

51* C U I = 0 E N 5 I T Y 

52* C U2 'MOMENTUM DENSITY _ 

53* C UU=VELOCITY 

59* C XaSPATlAL COORDINaTF MEASURED UPWARD 

55* C ZaS p ATIAL COORDINATE •' E A S \J » F D 0 0 ‘A N A R F.) 

S' C SOTsSHOCK our AM On time at the L U \ £ R BOUNDARY* 

57* READ (SUL- > f STOP 



b fl * „ 



READ ( 5*66 ) ( * T E'( f) UTTTT'n «T, 

$ 9 • 



I M | a I M A X - I 

60* 



JM2= IMAX-2 

6 1 * 



G AM = 5 * / 3 * — — — — - — 

6 2* 



GAM i «G AM- 1 . 

6 3,_ 



G AM3 * 3 • -G AM 

66 * 



R$s*695E+6 __ _ __ . 

65* 



HP=7 • E5 

66* 

♦ 


HC = 3.E^ 

67* 



G S ° • 2 7 R . 

69* 



R=.83E>2 

69* 



DTP*M8Q * 

7Q* 



SpTs6^Q* 

7 I * 



C 1 = G A M 1 / R 

7 2* 



GSR5= I G5*RS**2 

7 3* 



U1 0*»M. 1971 3E8 

7<i* 

C 



75 * 

c 


I N ITl A t ! Z A T ! 0 N 

76* 

c 



77* 



RX ( ( ) ■ P 5 ♦ HC 

7 fl * 



T = C* 

79* 



TP = Q. 

Bn* 



DO 20 1*1,1 

fl 1 * 



R X ( J ) = X 1 I ) * R S 

8 2* 



RXSQ f I ) =R X < 1 } * *2 

83* 


20 

■'CONTINUE 

BH* 



DO 77 I * 2 * I M A X 

8 S * 



X ( I ) " X ( I - I ) + 2«jGUl;« 

86* 



RX C 1 > =*X ( { J +RS 

07* 



RX50 < 1 > *RX ( i ) * *2 

88 , 



T E ( I ) 35 1 » 5 3 E 6 

B 9 * 



U|(1)*U10*EXP((GSrS/(R*TE(1)))*(K 

90* 


77 

CONT I NUE 

9 ] * 



DO 60 I * 1 » 1 MAX 

92* 



C ON « 2 * E - 1 2 

93* 



Ul ( ! )aCON*U] ( I ) 

9 H * 



UU< I >*0» 

95* 



U2( I )*Ut ( I )*UU< I ) 

96 * 



PH M M « + GSRS/RX < I ) 

97* 



U3(I>»U|C!1*<TE(I>/CI+tS*UU(t>**2| 

93* 



C ( I )«SORT ( GAM*R*Te ( 1 ) ) 

99 • 



El ( ! )«Ul C I )*TE U ) /C 1 

ICO* 



U 1 C ( 1 ) B u l ( I ) 

l C ! * 



U3C( 1 ) *U3( l > 

\ C 2 * 



tec m «te < n 

103* 


6C 

CONTINUE 

I 0 9 • 



DO J 1 1 1=2, (MS 

\ C5* 



H X ( I )=GSRS/(R*TE( [ I )•( (X( I + l > ~ X ( I ) 

106* 



CX(I>= <HX< J J **5* < EXP(-HX ( I ) )-CXP( 

f 0 7 • 


S I 1 

CONTI NUE 

108* 



WRITE ( 6 . I 0 ) T 

i 09* 



„ WR I T E ( 6 » 105 ) ... 

1 10* 



WRITE ( 6 » 1 I y ) ( X ( I ) , I ® 1 , I MAX ) 

111* 



WRITE ( 6 , l 05 J 

117* • 



WRITE (6, I 1- } (Ul ( I) » ! = ! » I MAX ) 

i • 



WRITE ( 6 i l 0 5 ) 

i 1 ; i * 



WRITE (6,110) (UU(I),!*J,IMAX) 

1 1 5* 



WRITE (6*105) 


/ R X < I ) » t o / R X ( j ) ) ) 


) /Rxsat i ) ) 

HX ( I ) ) > ) / < X < I + 1 ) - X U ) ) 


Q 



116 * 
i 17* 

1 \ 8 • 

I l 9* 

1 20 * 

1 2 l .lSr 

122* 

1 23* 

I 24* 

125* 

! 26 * 

1 27* 

I 2 8 * 

129* 

1 3,0* 

1 3 * * 

1 3 2* 

J 33* 

I 34* 

13 5* 

1 36* 
l 37 * 

1 3 8 * 

1 39 * 

ho* 

I 4 1 * 

* u 2* 

14 3* 

14 4 4 
» 4 5 * 
146* 
147* 

14 8* 

149* 

I SO* 

t 5 1 * 

i 52* 
i 53* 

i 5.4 * 

I 55 * 

1 56* 

i 5 7* 

! 58 * 

1 59* 

i 60 * 

1 6 1 * 

! 6 2 * 
i 6 3 * 

i 64* 
65* 

• 66 *.. 

1 67 * 

6 8* 

. 69 • 

7 0 
7 1 • 
72* 
73* 


WRITE (6,110) CU3 ( 1 ) % I *1 ♦ ! MAX ) 

WRITE (6, ICS) 

WRITE < 6 , II c ) ( T E ( I ) » I * 1 » I MJUO 

WRITE (6*105) 

WRITE ( 6 , | 1 y ) ( PH I ( I ) , ! * l | I M A X I 

WRITE ( 6 , I 2 £ 1 ' . . . 

c 

C BACKGROUND CORONA RADIATION rs DETERMINED -By TUC«ER»5 EQNo(>*«o 

c radiation effects are introduced by including radiation eqn_ 
c obtained by fitting of tucker *s results. 
c 

DO 62 I 3 1*IMAX , v 

I F ( T E ( I ) • G T ♦ 5 • u E 4 ) GO TO 11 . ' 

C2« 1 * jE|Q 
03=3*59 

""'GO TO 7 9 " "" - 

11 IF ( T E ( I ) • G T * 2 • 5 E 5 ) GO TO 12 

C 2* 3 » C'E 2 6 _ : ; 

C 3 = 0 * 

GO TO 79 

12 IF { T E ( I ) • G T . 7. tJ E 6 )_ GO T 0 1 3 _ _ _ _ 

C 2 = 1 ♦ LE 3 2 
• C3 = - 1 .172 
GO TO 79 

13 C 2= ! ♦ uE 2 3 ~ 

C3 *0*288 

79 QRC( I ) = C2* ( U I { I I • ♦ 2)* ( TE ( I ) * *C3 ) _ ' , 

62 CONTINUE 

c 

c _ 

c PLOTTING OF THE INITIAL VALUES' OF THE V A R I ArLES , * . « 

C 

CAL L. PL 0 T | ii. # _U L.UU.JJJ ; ^ 

TP=TP*OTP 
Ul I = U I < 1 ) 

„. U 2 I = U 2 ( 1 )_ __ 

UUfaUU(I) 

U 3 I * U 3 < 1 ) 

TE I » T E U i ■ „ 

C 

C 

C . ... upward, prop agat i ng shock... s i mul.at i on jn term s o r dens ity » 

C ENERGY , TEMPERATURE gradient..... 

c reflected shock m a ch numbers are assumed to be knck-jn 

c at. the lower boundary • then the.. Jump condi Tyo ns. are. pound by. 

c rank ine-hugoniot RELATIONS, 

c 

UIC = UJ. ( i 1*2. 

UUCi = 0 . 

U 2 D = U i 0 * U U 0 

... ... TEG* TE ( 1 ) 

U30»U l 0* ( TEtj/C l ♦. S*UU0**2 ) 

C 

KW* \ ' ...... 

U l ( 1 )*li!u 

uu< n = u u o 

u2 ( i )*u2u 

U 3 ( 1 ) * U 3 w " « 



ft ( I )VTE-u 

c( i )=sort<gan*r«te( j > ) 
£1 ( t )*U! ( I >*TE( J )/C! 

7 C CONTINUE 


PREDICTOR STEP 


. C ALL DT ! ME <UU,C ^--r r .-~r-r-rr 

IF COT. LT# a. on GO TO 200 
IF (T.LT.3?. ) DT*1. 

_.SDD*SDT*3Gf — — 

IF (T.GT.SDT.AND»T.LT.SDD) 0 T “ l • 

CALL FLUX I (UI •U2,U3».TE*UU I C) 

CALL E ON 1 < U J , U2 » U3 , F 1 , F,2.».F 3 t * 2 *A3 j PI »C X » .HX , K } f 

CORRECTOR STEP 

C ALL FLU X 2 < U I A # u J e , U 2 A , u 2 B I u 3 A , U 3 B • U 1 H , u 2 H g U 3 H 9 T E H D u U , c » 

I UUA , UUR ,UUH t UI) 

CALL E0N2 (uf,U2,U3,Fl,F2,F3,K? t K 3 * F I A , F 2 A » F 3 A * KJj> f 1 H , E l» 

lF|B,F2B,F3B,K2H,K3H,OT,CX,HX*KlA,KlB,K2A»K2BiK3A,K3BJ 

t*t+ot 

•- REFLECTED SHOCK SIMULA T I Q N ~ 1 N "T £ R M~ Or ~ T H E “K N OWN M ACH 'NU M B E R 


CALL BOUND „ 

I F ( K W * G T * G ) GO TO 55 5 

-• GO TO 666 

555 CONTINUE 

WRITE (6*101 T 

WRITE (6*105) _ 

WRITE (6* ML) (UI ( I> » 1*1 9 I MAX) 

WRITE (6,105) 

WR I TE < A , I 1 li ) { UU ( I > ,J ■ J. * _ I M_ A X J 

WRITE (6,105) 

WRITE (6,11c.) ( U3 ( 1 > « 1 * 1 * I M A X ) 

_W R I T E (6,105) 

WRITE (6,1 IQ) (TElli,I a l*lMAX) 
WRITE (6,135) 

WR I TE < 6 , 1 lu ) ( E I ( 1 ) , I » 1 » 1 M A X J 

WRITE (6,12L> 

00 ROM 1= I » IMAX 

C U )=SQRT ( GAM* R . * T £ (_I ) ) 

MOM CONTINUE 
KVV®0 

66 6 ...CON T 1 NUE 

KW*KW+ 1 


.PLOTTING O.T ... TH £ P R R T U R a £D..y AP I ABLES o . 0 . 5 . 04 , 


IF tT.GT.TP> 60 To l&O 
GO TO 160 

15 C CONTINUE 

CALL PLOT ( | ,U I ,UU, TE ) 
f pa TR + DTP 

1 6 0 IF { T,C,T, TSTOP ) Go TO 2C0 
GO TO 70 


0 



232 * 

2 33* 

234* 

2 35* 

236 
2 3 7 
23 8 * 

239* !' 

2 40* 

241 * 

2 4 2* 

243* 

244* 

2 4 5* 

246* 

.2 4.7 *_ 

248* 

249* 

250* 

25 1 * 

252* 

253* 

2 54* 

255* 

256 * 

257* 

258* 

_.2.5_5L* 

260 * 

26* 

262* 

263* 

264* 

26 5 * 

26 6* 

267* 

2 6 8 * 

269* 

27C* 

. 271 * - - 

272* 

273* 


2 0 C CONTINUE 
C 

K» ) 

DO 1-9-0 1*1 9 I H A X 
X ( K 1 * X ( I 1 

K = K* 1 
C 

I 9q CONTINUE 
C 

c the following is to record the. datas on the tape**;* 

C ' ~ ~ 

DO 2 I G J* I » JM AX 

0 0 2ll 1 = 1* IMAX _ _ 

IF U.EQ*IMAX> X ( I ) * 1 * oE 3 7- 
IF ( I . E Q * I M A X ) UIP<J»I)«1*0E37 

' I F < I . FQ * J H A X * A NO., J , F 0 , . jM.A X > . X (. | 1 = ? *0E37 .■/ 

IF <1 *FQ. IMAX. AND, J.EO.JMAX ) U 1 p ( J , I) a 9 • OF 3 7 
WRITE 191 X<l)«UlP(J«I) 

.210 CONTINUE 

DO 30 y J * 1 # JMAX 
00 30J 1*1 * IMAX 

if.. < r • e q • r m a x i x ( 1 1 * i • ae. 3 7 

IF ( I * E 0 , I M A x J U U P ( J , I 1 = i * C E 3 7 

IF ( 1 *EQ* I M AX . AND . J*EQ , JMAX ) X(1>*9,GE37 

1 F ( ! * E 0 • I M A X • A N D • j , E Q * jM A X ) UUP ( J t \ ) * 9 4 OE 37 

WRITE (9) X{ n ,UUPC J,I ) 

300 CONTINUE 

j. DO 4 00 J* 1. § JMAX 1 . 

DO 4 OL 1 = 1 * IMAX 
IF (I«EO*IMAX) X(f)**l*OE37 

IF— C ! *E0. IMAX 1--TEPC J » 1 1*1 *0E37 __ 

IF < 1 *EQ* IMAX* AND. J. EQ, JMAX 1 X(I)#9,CE37 

IF ( I *EQ. IMAX* AND, J,EO, JMAX 1 T E P I J , I 1 * 9 * Q E 3 7 

WRITE (9) X( M # .TEp( J, L1 , 

400 CONTINUE 
STOP 

JLC_ F_0 R m A t _ ( E 1.5 , 4 ! * 

66 FORMAT (2FI5.4,E15*6l 
105 FORMAT < 1 H ) 

lie. FORMAT (IX, 9E 11* 4 ) 

I 2G FORMAT < 1 H I > 

END 


ND OF COMPILATION 


NO DIAGNOSTICS 



I * 

2 * 

3* 

M* 

'&* 

6 * 

7 * 

a* 

9* 

!£.*_. 

11 * 

1 2* 

1 3* 
1M * 

IS* 

L6*.. 

17* 

18* 

1 9*_ 
20 * 

2 1 * 

22 1 
23* 

2 M * 
25*__ 
26* 
27* 
2_8*™ 
29* 

3 0* ' 
.3.1 *_. 
32* 
33* 
3M7L. 
3S* 
36* 
37-*- 
36 • 
39* 
. 40 *. 

8 l * 
M2* 

M 3 * _ 
MM* 
MS * 
.4.6* 

M 7* 
M 6 * 

M 9 * 
SO* 
S 1 * 

52* 
S3* 
5 M * 
SB . 
S 6 • 
57* 


SUBROUTINE’ FLUX 1 ( U I » U 2 # U 3 . T E * U U * C ) ^ 

PARAMFTFR IMAXcMt 
REAL K1 ♦ K 2 » K 3 

COMMON/fiLOC \ t 1 M 1 8 1 M2 _ r- 

C0MH0N/RL0C7/GAM »GAM! t GaH 3 .GSRS ,Cl»R 

COMMON/ RL0C3/RXC tMAX)tX(IMAX), PHI UMAX) iRX SO UHA > 

COMMON/ BL 0 C M / F I < 1 MAX) * F 2 M M AX) . F 3 U M AX > ♦* 2 J J MA X ) »A1< 
i .kin max) 

COMMON/ BLOC 8 /TFC I 1 MaX ) ,U lC I IMAX ) 

COMHON/BLOC I 6/OT 

C0MM0N/RL0CM2/QRIIMAX) 

COMMON/BLOCM3/QPC < t MAX ) 

DIMENSION U 1 U M A X ) , U 2 C I M A X ) . U 3 ( I M a X I t J XJ J M A X } 

DIMENSION UUUMAX)*cUMAX) # 00(tMAX)|DM(jMAX) 

DIMENSION ud< I MAX ) * .X... 


I M A X ) 


C2 AND C 3 ARE CONSTANTS USED 
CO*CONVER5 ION FACTOR 


IN TUCKER’S RADIATION EON * 


C0= ! •q£>22 
AC * 3 « 


GO TO 12 


" RAD I A T i 0 N E F> E C T 5 A RE INTRODUCED BY TnCLUOImG RADUTlON EQN 

OBTAINED RY FITTING OF TUCKER, S RESULTS* 

D o 6 u. I « l » t M A X - — r — : r- — — 

IF ( TE U ) .GT*5. jEm > GO TO II 
C2* fl « OF l u 

C 3 .» 3 * .5 5 — — — — — - - — 

GO TO 77 

It IF <TE < I ) .GT*2.SES ) 

0.2*3 • (IE 2 6 ..... — 

C 3 * 0 # 

GO TO 77 

L2- ! F I T E l I ) • GI. 7 .0,161 

C2»t«UF32 
C3--I • 172 

GO TO .7 7 — ~ — 

I 3 C2« I * 2 3 

C 3-0 « 288 

77 QRI I)*C2MU 1 1 11* « 2 1 * .1 TE-U J_*.*C3J*C 0 
60 CONTINUE 

DO 89 I ■ 2 i I M 1 


_G0_T..Q_1.3. 


80 


DO 89 I * 2 i I M l 

QRC I ) s-OR< !)♦.&*< * s*_( OR U * l >.♦9* U - 1 IX* Q R C LI .).« C.01 

CONT I NUE 

ORI I )»«5*I-QRI I )**S*IQR(2)^QRU ) ) ^ QR C I I ) 0 C 0 ) * 
QRUMaX)*.BU-QR(IMaX)*,5*JQRUM A_X1* GLR_U M 1 } J ± Q R C i-LM A 


66 


3 C 
2i 


00 66 I “2 t I M I 

U 0 ( I 1 * U U ( I ) - UU.t I -1 ) 

CONT t NUE 

DO 20 I r 2 * I M 1 

IF,. IUOU ) • GT *L * ) ..GO ..TO .30. 

GO TO 2 U 

Ut> < I )■ -0* 

CONTINUE 

U D 1 1 > a 0 * 

UD ( I M A X ) *0 • 



c 


$9* 

60* DO 50 ! °2 d I M S 

6 1 * ..... . DM<I )=UU( | • { )/C(! _) - uu ( I >/C( n ... — ~ 

( ;2 « I F ( DM < 1 ) • GT . 1 » a > GO To 80 

6 3 QO ( I )«CN 

. .6 H _i GO. TO . 5 0 _ - - 

65* 80 CONTINUE 

66* SM**OH<n 

6 7_* P A I s ( 2 • • G A M / ( G A M * J . ) ) • ( 5 M * * 2 > ~ ( 6 A M I / I GAM* I ♦ ) j 

68*~ SI G »( G A M*| * )*<SM**2> /(GAHl*SM**2*2« ) f 

69* ZETA B 5QRT<PAI/SIG) 

70* 00< I >=U1 (!)*<C(I>**2)*(«5*(PAT-l*)*(l«/StG*lt) 

~ 7 I * T-T GAM/GAM n • ( ZETA**2* I . I r/ (OT*GAM I " ~ 7 ~ 

72*. Q D < I ) a A B S ( Q D < I) ) ■ - <> . H 

_ 73* 5 C _ C ON T I NUE ..... ; 

7 8* QD(l)*0« 

75* 00 C I Ma X ) *0 « 

76 ♦ C 

77* DO 1C 1*1, IMAX 

78* F!(I)«U2U> 

79 * ^ V? t T )xGAMl*li3( 1 )*n»5*GAM3*U2( I I * * 2 / U I ( I ) : ' 

8 0* * " " 1 * • 5* (AC * *2 ) * U 1 ( I ) * ( UD I I ) * * 2 ) 

81* F3IMMU2in/U|{I)]MGAM*U3(I>-o5*6A^l*u2(l)**2/UlU) 

_ 82< L 1 * * 5 * ( A C * * 2 ) • U l ( I > ♦ ( U 0 ( I ) * * 2 ) ) 

83* K l < 1 ) a *2 o *U2 { I ) /P?X ( I > 

8 8* : - PC 2 t I ) a »U S < 1 )«GSRS/RXSQ( I ) " V: ^ V " - * ; ^ : 

6 5* 1 -2« oU2 M ) **2/ ( RX t H *U l < 1 > > \ v -P 

8 6* K 3(1 ) « - U 2 ( t ) * G SRS/r X $ 0 ( I ) 

87* 2-2o*U2(I)*(GAM«U3(I)-GAMl**S*U2(I)o*2/Ul(n 

8 8 3 * « 5* (AC * * 2 ) * U I ( I > * ( UD ( J ) * » 2 ) ) / < R X ( ]) * U l (JJJ 

8 9* *4*00(11 * 

90* IQ CONTINUE 

9 1* RETURN ■ ., ; P ' ■ • •' "" - 

92* ENO 

VCT OF “COM Pit A riON': “N0""D1AGN0S T ICST ~ ; --=-~r T -^r~ 


1 * 


SUBROUTINE eQNI ( U 1 . U2 f U3 , F I , F 2 , F3 , K 2 , K 3 * D T „ C X „ H X , K 1 > 


... ^ 

3 * 
u * 

5* 

6 * 

7 * 

8 * 
9* 

... 10 *.. 

1 I * 
12 * 

I 3 * 
18 * 

1 


parameter i m a x » <4 u v. v ; y: 

real K i 9 K 2 » K 3 * • : 

_C OMMpN/PLOC 1 / IM | * I M2 _ ' V / 

COMMON/OLOC2/GAMsGAM1 ,GA«3»G5RS 0 C1 8 R 
COMMO.N/BLOC3/RX ( IMAX) # X(IMaX) 9 PHI(IMAX),RXSq(1MaX) 

c 0 M M 0 N / R L 0 C | 1 / U * A ( I M A X ) » U 2 A I I M A X > 9 U 3 A ( I M A X ) , U I B ( I M A XJ * 

JU2B( IMAX ) * U 3 B ( IMAX) 

COMMON/BLOC I 2/ F U< IMAX ) *F2A c IMAX ) t F3A ( |MAX ) siFS BUMAX ) * 

! F 2 B ( I MAX ) , F 3 B J I MAX) ^ 

COMMON/BLOC1 3 / U 1 H < IMAX) 0 U2H( IMAx> 0 UUM( I M A X ) , T E H ( IMAX) 
COMMON/BLOC 1 8/K I H ( I MAX ) 9 K2 h ( IMAX ) 0 K3H ( | MAX ) 

COMMON/BLOC 1 5/U3H ( I MAX ) »UUA ( I MAX ) } UUB ( I MAX ) « 

DIMENSION HX(1MAX),CX(IMAX5 

DIMENSION Fl(lMAX) 9 F2(IMAX),F3(lMAX> 0 K2(IMAx) f K3(!MAX) 5 KltiMAX) 
DIMENSION U I (J max ) » U2 ( I MAX) 9 U3 ( I MAX ) 



IB* 

19 * 

2 G * " 
2l *w 

22 * 

2 3*' 
2 H* 

2 5 * 

2 6 * 
27* 

2 8 * 

29* 

30 ** 

3 i*_ 

32* 


00 ' 5 * ’ 1*2 , I M I 


uiAU) s . 5 *(ui(i*i>*uiUH'< 0 T/(K(i*n-xcm)*{n<r*i>- 

|,5*DT*IK{(I + l UK|(!) I — — -i— 

U 2 A ( I J«.S*(U2( f* l)*U2( I I )-t07/(X( !♦! )-X{II M*(F2C Wl )- 
t ( • 5 ♦. 0 T I * ( K 2 ( | + I . I * K 2 ( I I I _ , . 

U3A< 1 1 ®. 5 * r tn < i + i > ♦ u 3 * n > -J d t / c x < L 4 i .L* ^ i L!J_L* L? 3 !i *lL" 

I ( ,5 *DT ) * ( K3 { I + 1 > +K3 ( M > 

U|B(I) = .bMUl(I-U+U](l|) + <OT/IX(I)-Ml-n)|MF|lM) a 


F H I > ) * 

>777777 

F3 j m * 

•Fi i n ) + 


| • 5 * D T * I K 1 ( I - 1 I + K l { I > A _ - 

U 2 B ( I )«.5*(U2<f-l )*U2n))*(DT/(XU)-X(l 
| ( *5*PT ) * ( K2 ( I - 1 } +K2 ( I ) I 
U3B ( 1 ) * . 5 * l U 3J I - M * U 3 ( Ij > ♦ f p T / ( X t I > - X_< .X 


i ) ) ) * C F 2 < I - I > - F 2 U ) ) ♦ 
1 n ) •< F 2 ± I-l r»F3 ( 1 I ) » 


1 « .5*0T I * < K3 ( I - I ) + K3 ( 1 ) ) 

uiHt I >*•&* <ui u* n *un i- n x« i*i >*x< i -i > r-> * cr .1 u + n 

l-Fl (1-1 ) I + _ : 

2DT*K l < H 


3 3* 

3 q * ■ 

35* 

36* 

3 7*: 

38* 

39* 

BC*_ 

HI • 

M2* 

H 3*_ 
MM* 


2 + DT*CX < t ) *F i C I) 

3 ♦ • s * u i < r 3 * n .-.s* <expi-hx.( n > *exp (hx( i i j > > 

U 2 HU )«•>* fU2U* i I *U2( f-l > >-<DT / fX( I* ll*X f I-l ))) *(F2< 1*1 ) 
I-F2 I 1”1 ) ) * ; 

20 T * k ? ( T > .... _ - : ■■ ' ■ - 

2 + 0 T*CX C M *F2 U ) 

3 + * 5 * U 2 l 1 )* ( 1 *-. 5 * (EXP(-HX( 1) )+FXP(HX( I > ) ) ) 

U 3 H ( T ) = « 5 * ( U 3 U > I ) * U3 < I - l ) > - ( 0T / C X C I ♦ I ) * X ( 1 «> I ) ) ) » C F 3 < I * 1 I 

l-F3tM)> + 

2 D T * K 3 \ I ) • 7 • .., 

2 *DT*CX< I ) *F 3 ( I ) ...... ~ 

3+.5*u3e!>Mj.-.5*<exp(-HX<jn*EXP<HX(inn 


H 5 ' UUA ( n=U2A ( n /Ul a ( I ) 

M 6 U U S ( I ) = IJ ? 8 ( J ) / U I B ( I ) 

>7* ’ UUH ( I } *U2M ( I ) /U I H ( I I 

US* TEHU )=c I * (U3H( I ) / U 1 H < I >-* 5 *UUHf I|«*2) 

H 9 * • 50 CO NT 1 NUE i: ....; 

SO* ~ " RETURN 

5 I * END 


io op c6mpilation: no diagnostics# 



!♦ SUBROUTINE f L U X 2 ( U 1 A , u 1 B 0 U2 A , U2S 9 U3 A e U 3 8 s U I N 9 U2H s u 3 H , T EH , UU » C B 


2* _ t UUA , UUB i UUH , IJ U 

~ 3* PARAMETER iMAXsqy 

H * REAL K!H»K7h # K3H 

s» _ „REAL K 1 A , K I B , K 2 A » K 2B , K 3 A j K 3 B_._ 

6* ‘ COMMON/ BLOC I / l M ! # ! M2 

7 • COMM0N/RLOC2/6AM»r,AM» ,GAH3,GSRS,CJ P R 

8* - COMHON/0LOC3/RX ( IMAXI ,XUMAX ) ,PHI ( IMAX) »RX5Q( IMAX) 

9* " COMMON/RLOCS/TEC. c ! max ) , U l C ( I max ) 

jr>* COM MON / BLOC 1 2/ FlA{IMAX).F2AtlMAX),F3AUMAX) 9 F|BtlMAX)*, 

II . !F2B< I MAX ) t F 3 B t 1 M A X. I , . 


12 * 

1 3 * 

n* 

! 5 * 
16* 
17* 
18* 
19* 
20 * 
21 * 
22 * 
23f_ 
23 * 
25 * 
26 *_ 

2 7 * 
28 * 
29* 
30* 
31 * 
32* 


C0MM0N/BL0C I *4 /K 1 H ( IMAX)»K2H (I MAX) * K3H ( 
C0MM0N/BL0C25/KI A ( IMAX) «K I B ( I MAX ) *K2 A ( 

, K 3 8 ( l M A X ) . — 

COMMON /BL0C16/DT 
C0MM0N/BL0CM2/0B C 1 M AX ) 

DIMENSION U 1 A ( I M A X ) , U 2 A ( I M A X )_» U 3 A ( I M A X 
U 3B( t M A X 1 | U 1 H ( IMAX) » U 2 H < 1 M A X ) , U 3H < I M A X 
dimension UU ( I max ) ,UD ( 1 MAX ) 

c C 1 M A X ) , 0 M n M A X I , Q OJ I M AX) 

UD AC ! M A X > i M DBC IMA X»#UUAUMAX>,UUB(lMAXl 
UUH< 1 MAX ) 

U I ( I M A X ) ___ . 


IMAX) 

IMAX > ,K2B( IMAX ) »K3A< IMAX) 


) t U 1 0 
) i TEH 


( I MAX 

V J M AX 


) ,U2B { I MAX ) „ 

> " ’ 


66 


D I MENS 1 ON 
DIMENSION 

D I ME NS 1 ON 
dimension 
AC* 3 • 

DO 66 I * 2 * I M 1 

UD ( I )«UU C I )-UU{ 1 -1 ) 

U D A C I ) * C U U t 1 4 1 > • U U C n ) / 2 * 
UOB C U* CUUC ] )-UU< 1-1 ) ) /2* 

CONTINUE 

DO 20 l ** 2 • I M 1 

IF CUD C I 1 • G T • C * > GO TO 3 0 

6 0 TO 2 c 


33* 

30 

UD C I ) * 0 • 




33* 

2C 

CONT 1 NUE 




35* 


00. 21 1*2, IM) „ . 







36* 


IF CUDACn.GT.O* ) 

GO 

TO 

3 1 

37* 


GO TO 21 




3 3* 

3 1 

UOA ( I ) "U • .. . . 




39 * 

2 l 

CONTINUE 




3fl« 


DO 22 I *2 , I M 1 




3 f * 


IF CUDBC n ,GT.C. ) 

GO 

TO 

32 

3 2* 


GO TO 22 




M3* 

32 

U 0 B C I )*v* 




*4*4* 

22 

CONTINUE 



. — 

- 2 . 


35* 
*46* 
3 7 * 

Vn* 

39* 


DO 50 I* 
DM ( I ) * U U 

IF C DM C I 
OD C I ) »U • 


2 ♦ I M I 
C I - 1 ) / C l 1 

) « G T • 1 • 0 ) 


) -UU ( I > / C (J ) 
GO TO* 8 0 


.5.0* : 


GO TO 5u. — 

51 * 

80 

CONTI NUE 

52* 


5 M » D M ( l > 

53* 


PA 1 ■ C 2 . * GAM/ ( GAM* 1 , 

53* 


'SIG»<GAM*I« ) * < SM* * : 

55* 


2ETA»S0RTCPA1/SIG) 

56* 


Q pC I 1 aU| C 1 ) * CCC U - 

57* 


| - ( G A M / G A M l > * C 2 E T A « 

58* 


*QD( 1 ) * A B S C Q D I I ) ) 

59* 

5C 

CONTINUE 

60* C 



61 * 


DO 13 1 * 2 a T M 1 

62* 


F 1 A< I )=U2AC J 1 .. 

63* 


F 2 A C 1 ) “ G A M 1 o U 3 A ( I ) 

63* 


14.5*<AC**2)*IJ1A(I) 

65* . . 


F3A( I ) * ( U 2 A ( I ) / IJ 1 A 

66 ‘ 


l**5*CAC. « * 2 ) • U 1 A ( I ) 

6 7 * — 


FiB(!)*U28Cl» 

68* 

. • — . 

F2BC I ) « 6 A M 1 * U 38 ( 1) 


GAM* 1 * ) ) 


) » C 1 « /S 1 G+.l o ) 


5 * G A M .l * U 2 AC I ) * • 2 / U 1 A (.1) 



£9* ? ' I v. 5* f AC • *2 >"*U 1'B (T ) M UOR ( IT* * 21 


^ 0 * ’ F3Btn«(U2B( 1 >/UlB ( I ) > »<GAM«U3B( 1)-*5*GAM1 *u 2B< 1) **2/UlBU F 

7.1 * _ ] ♦ . 5 ♦ t AC**2)*U)B( U * CUDS til **2). ) 6 l 

72* K|H< I )»-2#*U2-H| ri/RX 1 II 

73* K2Hin»*UIH(n«(;SRS/RXS0<n 

71* 1 - 2 • « '1 2 H ( I >*•?/( R X < \ ) * U 1 H ( l± 

75* f(3HC I)*-U2H( I > *GSRS/RXSQ U ) , 

76* 2-2**U2M( l>*(GAM*U3H( I J ~G A M 1 * . 5 * U2H ( 1 >**2/UlH( II 

7 7 * 3* . S* t AC* * 2 )-* U IH ( UjM UD ( I ) R X 1 U .*U iHXlXl 

7 8 * <♦ ♦ Q D { I ) 

79* . K I A( I I*-2*«U2A( 1 ) /RX < I ) 

30*. K 1 B C I )*-2««U2B(.l > /RX.CI \ . „ 7 

3 1 * K2A { 1 M-U j A { I ) *GSRS/RXS0( I) 

32* !-2**U2A{ 1 )**2/(RX( I )*U]M I ) ) . ♦ 

33~* K. 2 B ( n a - U 1 B M ) * G s RS / R X S Q f 1 ) . li 

3 *1 * l-2**1J28(l)**2/(RX(T)*Ui8(l>) 

35* K3*< I )«-U2A(( )*GSRS/RXSQ< !) 

9.6 * 2 *• 2 • • U-Z A C I | * C G A M * U 3 A U ) - G A.M !..*. • 5 * U 2 A t I ) * *2/ U J A (JL> 

37* 4 3+«5*(AC**21«UtA< 1 1)»*2) )/<RX( I)*UlA(I) ) 

36* 4 K3B( n *GSRS/RXSQ( n 

9.9* 1 2 « 2 * * U 2 8 t n * I G A M * U 3 B 1 1) - G A M 1 ♦ .5 * u 2 a (I J A* 2/ UJL a.l L) 1 

90* 3*.5«<AC**2)*UIB(l)*(UD(I>**2n/(RX(I)*UjB(n) 

91* 10 CONTINUE 

92* R ETURN _ 

93* . END 

OF ~cT OMP 1 l a tTon: no dT ag NOS t ics « ~ : ' 




1* SUBROUTINE EQN2<U| ,U2.U3»Fl f F2,F3,K2,K3 f FlA,F2A t F3A t K| t KIH*EI, 

?* _____ lF|B t r2B f F3e t K2H > K3H t DT,CX,HX i K|A t KlB l K2AjK2B > K3A,K3 8l 

3* parameter i m a x » r j _ 

4* REAL K ! » K 2 » K3 

\» _ REAL KlH,K2H f K3H 

S* " ~ REAL K1 A . K S B « K 2 At K 2R, K 3 A *K3B ' ’ 

7 * COHMON/BLOC I / I to 1 « I M2 

3* COMMON/BLOC2/GAM»GAMi 9 GAM 3 , GSRS ,CJ S R 

7 * “* " " " C0MM0N/BL0C3/RX(lMAX),XtI M A X > , P H I" ( IMA X ) , R X S Q < I*M A X f “ 

3* CqMM0N/BL0C6/U3vi 9 TEil »U2D » uuo » U I o 

L* COM MO N / B L 0 C 7 / U I l ,tJ-2| .UUJ » U 3 I ? _T E I ■ . /. = _■ ,-i 

i * C0MM0N/BL0C8/TECI IMaX ) ,UJC( JMAX ) 

3* COMMON/BLOC2L/UUf IMaX ) , TE ( IMAX ) 

^* COMMON /BL0C33/U3C(IMAX)|UUC(IMAX) 

a* ’ DIMENSION HX ( I MAX ) f CX ( IMAX ) ~ ‘ 

i* DIMENSION F|{IMAX) # F2(lMAXl 9 F3(!MAX),K2C!MAX>*K3nMAX| t KlClMAX) 

7«_.; DIMENSION F1A.UMA X ) t F2A C, I.M A XJ * F3 A ( J M A X J 6 F 1 B ( IMAX) ,F2BJ ]MAX$ # 



l R * 

| 9 * 

i a* %, 

7 1 f 

22* 

_ 2 3 * 

7 M * 

25* 

2$ * 

l F38 t I M A X >, K 2H C 1 M A x ) , KTH ( IMA X ) t K J H M M A X ) 

DIMENSION Ul ( I M A X ) ,U2C I M A X 1 » U 3 ( I M a X ) 

DIMENSION K I A C I M A x 1 , K 1 B ( I MAXi_j.K2 A JLJ MAX.) # K2BjJ MA X ) 

I ,K3A ( IMAX ) ,K3B C IMAX ) 

0 1 MENS ! ON E I ( IMA X ) 

DIMENSION E I T I I MAX ) ,E IE ( IMAX ) 

oo mb 1 = 2* t m i ' -rV ; '* .•••,?■• i' ' 

mb r i r ( n =e i { i ) 

fin 50 I ■ 2 f ! M i ■ ' ... :■ - : ' u ;; -v-V./VI. 

27* ; 

28* 

2 , 

3 C * 

3 1 * 

3 ? * 

Ul(I)cUlM>-(DT/IX(I* , n-XCI-I)))*C«5*IF|II*l)-Flt!-n) 
I ♦ ( F 1 A ( I) ~F 1 8 C I ) > ) 

2 +DT**S* (K I < I ) *K *H< ILL . . 

2* DT*CX C I) *F I ( I ) 

U 2 <!>*U 2 (!)-CDT/(X(I + n«-XII«*in)M*&*(F2lI + l>-F2tI-l)> 
1*(F2A{P-F?P(I)>) - 

33* 

■ 34* 

... 3 S.* __ 

3 A* 

37* 

.38* _ 

39* 

MO* 

M ) * 

2*DT* • S* ( K2 ( I > +K2H ( I) > 

2* D T * C X f I ) *F 2 I 1 ) 

_ U3 « .1.) «U3 I I > - < D T / C X ( L*U»XII« F_3Jl*_lL-F3(I-lL) — 

1 + CF3A II >-F3B ( I ) ) ) 

2 + D T * • 5 * l K 3 C I ) ♦ K 3 H ( 1) ) 

2*DT ♦ C X Cl ) *F3 l I ) 

UU< I > *U2 < I ) /Ul ( t ) 

TEC I 1 a C 1 * ( U3 C I )/Ul C I ) - • S* UU C I )**2) 
EIU)=U1!1)*TEU) /.Cl 

M2* 

50 CONTINUE 

M3* 

C ' . 

MM* 

DO M 7 I a 2 » IAl _ ■ ■ -Z-ZiLL 


MS* EIE(1WUT< I I-ElC n )/ElT< |) 

IF (E1F ( I ) .LT .C.SO ) GO TO m7 

_M7* El C I ) * Q • 5 0 * E I T (I ). 

MR* TEC l i*El < I i *C I/Ul U ) 

U3< I ) « t.l 1 C ! ) • ( T E ( I )/CI*»5*UU< 1 1**2) 

...50 * M 7 C 0 NX I N U E - 

5 J * C 

S2* RETURN 

5 3 , * - END _ 


[5 OF COMPILATION 


no DIAGNOSTICS, 



) * 

?.* 

3* 

M*_ 

5* 

6 *' 

7* 

R* 

9 * 

. 

1 * 

2 * 

3* 

M* 

5 * 

6 * 
?*“' 
B * 

!G * 

! 1 * 

> 2 .*~_ 

?3* 

> q * 

2-S.*_ 

26* 

27* 
? 8 *_ 
29* 
30 * 

3 l.*_. 
32* 
33* 
3M,*._ 
35* 
36 * 
.3Z*- 
38* 
39* 

m.g*.._ 

H 1 * 
M2* 
__S 3 *- 
MM* 
MB * 
H6 *_ 
M7 * 

M 8 * 
M9* 
SO* 


Subroutine bound _ 

PARAMETER ! m A X * MO 
COMNON/SLOC 1 / 1 Ml » !M2 

COMMON/ BL 0 C 2 / 0 A H » G A M 1 , 6 A M3.».G SR S t C 1 » R 
COMMON / BLOCS/U! ( l M A X > , U 2 * IMAX) » U 3 < IMAX) 
COMMON/BLOC 6 /U3C • TEG * U2Q i UUG # U 10 
COMMON/BLOC 7 / U 8 I * U 21 i U U I t U 3 1 j.TJEJ_ - 

COMMON/BLOC8/TEC<lMAX) # UJCttMAX) 

C0MM0N/BL0C2C/UU(lMAX),TE(|MAX) 

CO M M 0 N / B L o C 3 3 / U3C ( I M AX J_*.U 

common/blocsc/t 
C0MMON/BL0C3^/E!< i m ax j 


SDT»6C0» — — 

IP CT.GT.SDT) GO TO 10! 

um >»u!Q 

U2 < M sU2 C 2 ) 

• ■ ” uu< 1 i = u 2 1 n/uiu) 

TE( 1)=TEC« 

El Cl) = U 1 <1 ) * TEX! > /Cl 

U3( I >“U1 { 1 ) *(TEI 1 )/CI+#5*UU( l )**2) 

GO TO 851 

l.Q I JC 0 N T I N U E ... T-.- 

u 1 ( ( ) =U t 1 
T£C I >»TE! 

U2 l l ) *0* '• - — — : : — — ■■ 

uuc i ) =U2 n > /u i n > 

E 1 { 1 ) =UH I > *TE ( 3 ) /C » 

U3< i )*U1 ( I >* <T E (U / C 1 _*_t.5 a U ULLl* ® 2 1 


IS1 CONTINUE 
C 

C U P.PE R - BOU N D A R.Y—V A L UE._t_§. • — 

C 

IF C UU ( 1 M 1 ) , GT # 20 « > GO TO |S0 

U i « 1 M A X ) *U I Cl I MAX) r 

UU< IMAXI«UUC< IMAX) 

U2tlMAX)aUt(lMAX)«UU(IMAX) 

TXII M A X.) * T ECC IMAX i ~ ~ 

U3 t ! M AX ) «U3C ( I M AX ) 

E!flMAX)*UK!MAX)*TEUMAX)/Cl 

.GO... .TO 105 : — 

ISO U1 i 1MAX)=U1 < 1 M I ) 
U2<1MAX1=U2(|M|) 

U3 i 1 M A.X ) »U3 U M U — — - 

UU 1 1 M A X ) SUU { 1 Ml ) 

♦TE< 1MAX)*TE( I M I ) 

Eli IMAX ) *F It I. Ml > — 

105 CONTINUE 
C 

.RETURN — •— 1 - 

END 




D 0 r C 6 MRILAT 10 N 


NO DIAGNOSTICS#' 



’ " T" 

2^- 

3* 


SUBROUTINE "PLOT ~ H FL'AGiU 1, UU , tE 1 
PARAMETER IMAXMU 

PARAMFTFR JMAXaG 

, H* 


COMMON/8LOC9/UI P( jMAXVTmAX ) ,IJUP{ JMAX, IMAX) ,TEPt JMAX* IMAX ) 

5* 


DIMENSION ui ( IMAX) » U U < IMAX) , T E ( I M A X ) 

6* 


!F ( TFL AG. EG. 1 ) GO TO 20 

" ”7* 


J~ 1 

8* 


K* J 

9* 

20 

CONTINUE _ 

10* 


DO 10 I - 1 ♦ I M A X 

1 1 * 


U) P ( J,K )=u 1 ( M 

I 2* 


UUP ( J , K ) «UU ( p ,:;1 il — * ■ . . . . 

13* 


TEP C J,K ) -TE f I ) 

18* 


K*K+ I 

1 5 * 

1 0 

CONTINUE 

16* 


J = J+ 1 

17* 


K = 1 

1 8 * 


RETURN- ' ■ - - • - - Z$.4i r ■■■ "■ •> . <A ,v> 

19* 


END 


rn or comp i l at i on : "no o i a gnostics “ 



\ 


\ 

I 


" | * " 5 U B R O U T I N E D T l M E ( U U » C > 

2* PARAMETER IMAXs^C, 

3* . COMMON/8LOC I / I M \ , I M2 , 

*4* COMMON/RLOC2/OAM f GAMl * G A M 3 » GSRS , C 1 ,R 

5* COMMON/BLOC3/RX( IMAX) ,X( IMAX) ,PHl C IMAX) t RXSQ( IM.AX1 

6* COM MON/RIO C | 6 / D T ■ . , . ■ , 

7* DIMENSION UU(IMAX), C ( I M A X ) 

8* DIMENSION DOT UMAX) 

DO 1 3 I * * » I M l .. : : . 

10* DOT ( I )=,2&* ( X< I + Ij-X( I ))/< ABS(UU( 1) >+C{ I M 

I 1 * 10 cont inue 

1.2* . H D T * 0 D T ( 1 ) - ■ . ' ' . . ' , ^ •• •• »' 

1 3* 00 20 ! “ 2 * I M I 

| 8 * GOT sDDT ( I ) 

._JL5 * : IF ( HOT , GT , QpT ) G 0 ..T.Q 3 a 

! 6* , GO TO 20 

17* 3 C HOT = GDT 

_ 1 8* 2 Q, C 0 N T I N U t_ • - . 

|9* D T = H 0 T 

20 * RETURN 

21* _ . END 


io of compilation: no diagnostics. 



APPENDIX C 


Presentations 


During the period of performance of this contract, the following 

papers were presented. 

(1) "A Compressible MHD Model of the Development of a Sunspot," 

Annual Meeting of the Solar Physics Division, American Astro- 
nomical Society, Huntsville, Alabama, November 17-19, 1970, with 
M. Hagyard and Y. Nakagawa. 

(2) "H Flares: The Response of the Chromosphere to a Downward Shock 
Wave," Annual Meeting of the Solar Physics Division, American 
Astronomical Society, University of Maryland, April 4-6, 1972, 
with S. M. Han and Y. Nakagawa. 

(3) "Non-Linear Study of the Dynamical Behavior of Force-Free Magnetic 
Field," Annual Meeting of the Solar Physics Division, American 
Astronomical Society, University of Maryland, April 4-6, 1972, 
with M. Hagyard and Y. Nakagawa. 

(4) "Some Characteristics of Disturbed Solar Atmosphere," High Altitude 
Observatory Seminar (Invited), July 1973. 

(5) "Solar Atmosphere," AIAA (American Institute of Aero. & Astro. )- 
Alabama Section Space and Atmospheric Sciences Panel Meeting, 
(Invited), December 3, 1973. 



